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Abstract 

Recently developed series representations of the Boltzmann operator are used to obtain Quantum 
CN | Mechanical results for the matrix elements, (x\ exp(—(3H)\x'), of the imaginary time propagator. 



The calculations are done for two different potential surfaces: one of them is an Eckart Barrier and 
the other one is a double well potential surface. Numerical convergence of the series are investigated. 



Although the zeroth order term is sufficient at high temperatures, it does not lead to the correct 

O ' saddle point structure at low temperatures where the tunneling is important. Nevertheless the 

ce 

O ■ series converges rapidly even at low temperatures. Some of the double well calculations are also 

J>^, done with the bare potential (without Gaussian averaging). Some equations of motion related 

Qh| with bare potentials are also derived. The use of the bare potential results in faster integrations 

of equations of motion. Although, it causes lower accuracy in the zeroth order approximation, the 

> '. 
(j\^ . series show similar convergence properties both for Gaussian averaged calculations and the bare 

CT\ potential calculations. However, the series may not converge for bare potential calculations at low 

in ' 

-^ ■ temperatures because of the low accuracy of zeroth order approximation. Interestingly, it is found 

£\j . that the number of saddle points of (x\ exp(—/3H)\x') increases as the temperature is lowered. An 

explanation of observed structures at low temperatures remains as a challenge. Besides, it has 
K> ' implications for the quantum instanton theory of reaction rates at very low temperatures. 
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I. INTRODUCTION 

Semiclassical theories are useful for the treatments of large dimensional systems because of 
their favorable scaling with the dimension of the system. They allow one to study systems 
that are very hard to treat with quantum mechanics due to exponential scaling of the 
sizes of bases with increasing number of dimensions in quantum mechanics. Besides, they 
offer advantages over classical methods by allowing one to include quantum effects such as 
interference and tunneling which cannot be described with classical methods. 

Development of semiclassical propagators for studying dynamics of atomic and molecular 
systems dates back to pioneering works of Van Vleck [lj] and Gutzwiller |2|. Van Vleck 
Propagator is obtained as a semiclassical approximation to path integral and it is exact 
for quadratic Hamiltonians. However, it has some drawbacks due to some difficulties re- 
lated with its numerical implementation. Firstly, dynamics is formulated as a double ended 
boundary value problem, so that it is necessary to do nonlinear searches to find classical 
trajectories that satisfy the boundary conditions of the problem. Secondly, Van Vleck prop- 
agator has a prefactor that has singularities. These problems with the Van Vleck propagator 
is overcome by the development of initial value representation (IVR) methods. In the IVR, 
integration variable related with the end point of a coordinate is transformed to an initial 
momentum so that the dynamics problem is posed as an initial value problem in which the 
initial conditions of the problem are specified by the initial phase space points. Thus, it 
is not necessary to make nonlinear searches for finding classical paths connecting the end 
points. Besides, the IVR based propagators do not include prefactors that have singularities. 
In addition to its numerical advantages IVR also offers a more intuitive physical picture for 
dynamics. 
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History of IVR dates back to its use by Miller [3[ and Markus |4[ in studies of classical 
S-matrix calculations for collisions. Modern semiclassical propagators are based on the idea 
of using Gaussian wave packets suggested by Heller p, |6] which is later refined by Herman 
and Kluk [7[. Since the work of Heller, many IVR based propagators are developed for real 
time dynamics. Several reviews about different aspects of the subject can be found in the 
literature [8r|l6|. 

The success of the IVR based methods in real time dynamics also motivated the de- 
velopment of semiclassical methods for imaginary time dynamics. Several methods has 



been proposed for a semiclassical approximation of the imaginary time propagator 17H20]. 



Recently^ 



agator 



19 



rantsuzov et. al. developed semiclassical approximation to imaginary time prop- 



20|, named Time Evolving Gaussian Approximation (TEGA), which is based 
on an earlier method suggested by Hellsing et. al. 17j. Another approximation to imagi- 



nary time propagator is suggested by Pollak and Martin-Fierro [2l| . In this method, named 
PSTEGA (Phase Space Time Evolving Gaussian Approximation), coherent states are used. 
Since the Gaussians that are used in the TEGA method can be considered as coherent states 
with zero momentum, PSTEGA method can be considered as a generalization of the TEGA 
method. Although the propagator involves phase space integration instead of configuration 
space integration, it is possible to integrate equations of motion in an efficient way. As 
shown in the appendix, the momentum degrees of freedom can be integrated implicitly so 
that matrix elements of the equilibrium density matrix can be evaluated with an expansion 
in configuration space as in the TEGA method. Besides, the PSTEGA method provides a 
new way of evaluating time correlation functions. 

Although the IVR based semiclassical methods has been used successfully in many sys- 
tems, the approximations involved in these calculations remained uncontrolled such that 
there was no way to estimate the errors in these calculations. This problem has been over- 



come by the development of correction operator formalism by Pollak and co-workers [22H27|. 
They have shown that the exact quantum mechanical propagator can be expanded in a se- 
ries in which the zeroth order term is an approximation to the exact propagator. Then, the 
higher order terms are obtained from the zeroth order term in a recursive manner by using 
the correction operator. Therefore, it is possible obtain quantum mechanical results, at least 
in principle, starting with a semiclassical calculation. 

In this paper, some equations of motion related with the PSTEGA calculations, in which 
bare potentials are used, are derived. In addition to that an efficient way of solving the 
PSTEGA equations of motion is shown. Then the TEGA and the PSTEGA methods are 
applied to two different systems one with an Eckart Barrier and another with a double well 
potential surface. Numerical convergence properties of the series are investigated and the 
accuracy of the results is checked with quantum mechanical calculations. In the calculations 
both the Gaussian averaging and the bare potential are used. The results show that the 
use of bare potential leads to lower accuracy in the zeroth order approximation. Neverthe- 
less, calculations show similar convergence properties with the Gaussian averaged potential 



calculations. The results show interesting structures at low temperatures. 

In section [Til theory of the TEGA and PSTEGA methods are reviewed. Then, in section 
UTTl calculations and results are presented. The paper ends with discussions and conclusions 
in section IIVI 



elements of the Boltzmann operator 
Pollak and Martin-Fierro 



II. THEORY 

The theory of the TEGA method as a semiclassical approximation is first developed by 
Frantsuzov and Mandelshtam [19|, |20( based on an earlier method suggested by Hellsing 
et.al. [17|. Later, Shao and Pollak applied the Correction Operator formalism to the TEGA 
method and showed how it can be used to obtain quantum mechanical results for the matrix 

The theory of the PSTEGA method is developed by 
2l| . In the following subsections, first the theory of the correction 
operator formalism and the series expansion of the thermal propagator will be given in a 
general manner. Then, the details of TEGA and PSTEGA methods will be given. Equations 
of motion for the PSTEGA calculation in which the bare potential is used will be given for 
the first time. 

A. Preliminaries 

Consider an N dimensional system with the Hamiltonian 

H = ^ + V(q), (1) 

where q and p are N dimensional vectors of mass weighted coordinate and momenta respec- 
tively satisfying the usual commutation relation [qi,pj] = ihSij, and V(q) is the potential 
surface of the system. 

The thermal propagator K{r) = exp(—rH) is the solution of the imaginary time 
Schrodinger equation (or the Bloch equation), 

-§f ~ H) k ^) = °> (2) 

at imaginary time r with the initial condition K(0) = /, where I is the N dimensional 
identity matrix. 



If there exist a good approximation K (r) for the exact propagator K(t), then the cor- 
rection operator C{t) can be defined as follows: 

C{r)={-^--H^k Q {r). (3) 

The differential equation above can be inverted to an integral equation by realizing that the 

exact propagator is the solution of the homogeneous equation (Bloch equation). The formal 

solution is given by, 

K(r) = K (t) + / dr'K(r - t')C{t'). (4) 

Jo 

B. Series Representation Of The Thermal Propagator 

Although equation (j4j) provides a formal solution, it is not very useful in that form since 
the exact propagator appears on both sides of the equation. On the other hand, if K (t) 
is a good approximation to the exact propagator, K(t), then it makes sense to expand the 
exact solution in a series where Kq(t) is the leading order term of the series such that 

oo 

K(r) = J2Ur). (5) 

i=0 

By plugging the expansion above to equation (J3J) , and assuming that Kj ~ & , the following 
recursion relation is obtained for the higher order terms by equating the terms that are of 
the order of the same power of the correction operator: 

K l+1 (r) = r^Kiir - t')C(t'), i > 0. (6) 

Jo 

Thus, given an approximation, the exact thermal propagator can be obtained from that 
approximation recursively, by using a series in which the zeroth order term is the approxi- 
mation. 

C. Symmetric Form Of The Series Representation 

The exact thermal propagator is Hermitian: K(tY = K{r). However, both the TEGA 
and the PSTEGA approximations do not provide a hermitian representation for Kq{t). 
Therefore, the series expansion that contains the TEGA or the PSTEGA approximations 
as the zeroth order term cannot give a Hermitian representation of the exact propagator. 



This problem can be remedied as follows. In order to make every single term in the series 
expansion Hermitian, the equation below can be used: 

K(r) = K{t/2)K\t/2). (7) 

Since (K(t/2)K(t/2)^ = K{t)K{t/2)\ the use of this identity guarantees generation of a 
Hermitian representation of K(t) from any representation of if (t/2) regardless of whether 
that representation is Hermitian or not. Thus, by expanding all of the terms in equation (JTj) 
in a series as in equation flS]), the following series expansion is obtained for K(r) in terms of 
the terms in the series expansion of K (r/2) 

K(t) = J2k^(t), (8) 

3 

where 

k^(T) = f^K,(r/2)K^(T/2). (9) 

4=0 

It is clear that a representation of the term on the left hand side will be Hermitian regardless 
of whether the representations of the terms on the right hand side are Hermitian or not. 

D. Definitions Of Averaged Quantities 

Frantsuzov and Mandelshtam derived the equations of motion for the TEGA method 



variationally 20(. This method leads to the result that the potential and its derivatives 
should be Gaussian averaged. In the rest of the paper, the following notation is used for 
denoting the Gaussian averaging of a quantity h(q): 

/ i \ N / 2 i roo 

(h(q)) = (-) , \ dx 

W v/|det(G(r))| loo 

x exp(-(* " q{r)) T G{r)-\x - q(r)))h(x), (10) 

where q(r) is an N dimensional vector defining the center of the Gaussian and G(r) is an 
N x N dimensional positive definite matrix defining the width of the Gaussian. 

Shao and Pollak has shown that the equations of motion for the TEGA method can be 
derived by requiring that the method is exact for harmonic potentials 28| , and the same idea 



is also used in the development of the PSTEGA method 2l|. By using this idea of deriving 



equations of motion, Shao and Pollak suggested that the TEGA method can be generalized 
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for an arbitrary averaging function. Given a reasonable zeroth order approximation, the 
series representation will converge to the correct result. Consequently, Shao and Pollak 
suggested that the calculations can be done with the bare potential (without Gaussian 
averaging). Provided that this still gives a reasonable zeroth order approximation, one can 
obtain the exact propagator from that approximation via correction operator formalism. 
The use of the bare potential leads to faster integrations of equations of motion since the 
Gaussian averagings are not done. 

E. Equations Of Motion For The TEGA 

Matrix elements of the approximate solution of the Bloch equation in the TEGA is given 



by (l9|, [20] 



\ N/2 
(x\k (r)\q ) = ( — ] 



2nJ | det (G(r) ) | V2 
x exp (-!((* - q{r)) T G{T)- x {x - q(r))) + 7 (r)) . (11) 

In the equation above, the N dimensional vector q = q(r) defines the center of the Gaussian, 
the N x N dimensional positive definite matrix G = G{r) defines the width of the Gaussian, 
and the parameter 7 = 7(1") is a real scale factor. 

In order for this representation of the thermal propagator to satisfy the initial condition 
that K(0) = I, the following conditions should be imposed for small r: 

q(r ~ 0) = q , G(r ~ 0) = h 2 rl, 7 (r ~ 0) = -rV(q ). (12) 



If Gaussian averaging is used, the equations of motion for the three variables are 20l.l29|: 

±G(t) = -G(r)(VV r y( g (r)))G(r) + ^/, (13) 

±q(r) = -G(r)(VV(q(r))), (14) 

A 7(r) = d T r[(VV T y(q(r)))G(r)] - (V(q(r))}. (15) 

The use of the bare potential leads to the same equations of motion except that the potential 
and its derivatives are not Gaussian averaged. However, the coefficient —1/4 in the equation 



of motion of the variable j(t) becomes —1/2 



28 



When equation (Q is used the matrix elements of the each term in the series expansion 
of the Boltzmann operator can be calculated as follows: 



«=j 



i=0 



x\K^(r)\x') = J2 / dy{x\k l {T/2)\y){y\ky i {r/2)\x'). 



(16) 



F. Equations Of Motion For The PSTEGA 



In the PSTEGA method 2l|, the thermal propagator is represented in a coherent state 
basis whose coordinate state representation is given by 

_l/4 

(x\g{p,q,r)) 



det(G(r)) 



IX 



N 



exp (- l -{x - q (r)) T G-\T)(x - q(r)) + % -p T {r){x - q(r))\ . (17) 

In the equation above, N dimensional vectors q = q(r) and p = p{r) are the position and 
momentum vectors respectively and the N x N dimensional positive definite matrix G(t) is 
the width matrix. The matrix elements of the imaginary time propagator can be expanded 
in the coherent state basis as 

dp&q 



[x\ exp(— tH)\x') 



2n 



(x\ exp{-rH)\g(p, q, 0))(g{p, q, 0)\x'), 



(18) 



and the mixed matrix elements of the thermal propagator are approximated as 

(x\exp(-rH)\g(p,q,0)) ~ (x\k o (r)\g(p,q,0)) 

= f(p,Q,T)(x\g(p,q,T)). 



(19) 



If the potential is Gaussian averaged, equations for the variables, q, p, G, and /(r) are 
;iven by 21] 

dq{r) 



Or 

dp{j) 

dr 

9G(r) 

dr 



-G(r)(VV(q(r))), q(0) = q , 

-ft 2 G(T)-V(r), p(0)=po, 

-G(r)(\/V T V(q(r)))G(r) + h 2 I, 



f(r) = ex P (-£dr' 



+ ^Tr[G(r'y 1 ] 



P 1 (t')p(t') + (V(q(r'))) 



n dq(r') 



h P[T) - dr> 



(20) 
(21) 
(22) 

(23) 



As the way Pollak and Martin-Fierro have done 2l|, the equations of motion for the bare 
potential calculations can be derived by requiring that the method is exact for quadratic 
potentials. In this case /(r) is given by 

"1 



/0 



exp 



dr' 



p T (r')p(r') + V{q{r')) 
.„ dq(r') 



+ ^Tr [G (r>r]-- h P^r,- dT/ 



TT[VV T V(q(r'))G(r')] 



(24) 



Other equations will be the same except that the potential and its derivatives are not Gaus- 
sian averaged. In the PSTEGA approximation, initial width of the Gaussians is arbitrary. 

In calculations, it is not necessary to make a propagation in phase space, since the 
equation of motion for p(r), equation (I2~TT) . can be integrated implicitly. For a detailed 
explanation of how to integrate equations of motion in an efficient way see the appendix. 



G. Matrix Elements Of The Correction Operator 



Matrix elements of the correction operator are given by 21 



(x\C(T)\q(r)) = -{V &nh {x,q,T)){x\K Q {r)\q{r))., 



(25) 



where (V r an h(aJ, q, r)) is the anharmonic remainder of the potential when it is expanded 
around q(r). If Gaussian averaging is used, expansion of the potential surface is given by 



V(x) = (V(q(r))) + ^(V T V(q(r)))(x-q(r)) 

+ (x-q(r)) T (VV(q(r)))) 

+ l -{x-q{r)) T (VS/ T V{q{r))){x-q{r)) 

--\iyV T V{q{r)))G{r)] 



+{V anh (x,q,r)). 



(26) 



If the bare potential is used in the calculations; then, the expansion of the potential surface 
is given by 

V(x) = V{ q {r))+ l -{V T V{q{r)){x- q {r)) 
+( x - q ( r )fVV(q(r))) 
+ 1 -(x- q (r)rVV T V(q(r))(x- q (r)) 
+V anh (x,q,r). (27) 

In this case, matrix elements of the correction operator are given by 

(x\C(r)\q(r)) = -V^ h (x,q,r)(x\K (r)\q(r)). (28) 

H. Quantum Mechanical Calculations 

In order to make a comparison of TEGA and PSTEGA calculations with a direct quantum 
mechanical calculation, quantum mechanical calculations are also performed. In order to 
calculate the density matrix elements, first the Hamiltonian matrix is diagonalized; then, 
the thermal propagator is expanded in the basis of the eigenstates of the Hamiltonian. If 
4> n {x) are the eigenstates of the Hamiltonian which are obtained as a result of diagonalization 
calculation; then, the matrix elements of the thermal propagator can be evaluated as follows: 



•'"' 



'\exp(-PH)\x) = ^(a:'|exp(-^)|0 n )(0 n |x) 

n 

= 5> n (x')exp(-/3£ n )0;(^)- (29) 

n 

This calculation, in the case of Eckart Barrier, duplicates the calculations of Miller et. al. 

Q. 

III. CALCULATIONS AND RESULTS 

The calculations are done for two different potential surfaces. One of them is an Eckart 
Barrier and the other one is a double well potential surface. The details of calculations and 
their results are given in the following subsections. 
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A. Eckart Barrier 



The configuration matrix elements of the Eckart potential were previously studied by 



Miller et. al. 30|. The asymmetric Eckart barrier potential has the form 



v(x)= v.d-«) ^w+f; po) 

l + exp(-2ax) 4cosh 2 (ax) 

ren a = 1. 



30[, that is 



where a is the asymmetry parameter. The potential surface is symmetric w 

In the present study, the same parameters that were used by Miller et. al. 

Vq = 0.016a.u., a = 1. 3624a. u. and m = 1061a.u. are used. Some computations for an 

asymmetric barrier with a = 1.25 are also performed. 

Since this is a one dimensional system, it is convenient to work with matrices of the 

zeroth order Boltzmann operator, equation (II ip . and the correction operator, equation fl25|) . 

in the configuration space. Then, the final results are obtained by matrix multiplications and 

;ime integrations. The latter were all performed using the third order Simpson integrator 



3l| . An evenly spaced grid is taken for the coordinate y in a finite symmetric range, and 
a Gaussian form is defined around each grid point that satisfies the initial conditions given 
in equation ( JT2"|) . For the calculations presented in this paper 200 evenly spaced grid points 
in the range (—8, 8) is sufficient for converging the configuration space matrix elements of 
the Boltzmann operator in the range (—6, 6). For calculating time dependent averages as in 
equation (TTOl . Gauss- Hermite quadrature is used. 

To obtain the configuration matrix elements of the Boltzmann operator exp(—/3H), the 
equations of motion, equations (fl3"|) - ([T5"|) . were integrated up to the half time h/3/2 using 



the adaptive step size Cash-Karp Runge-Kutta method [3JJ. The matrices of the TEG A 
propagator and the correction operator are calculated and stored at every time step. Then, 
the higher order terms in the series corresponding to half time are calculated recursively 
by using equation ([6]). Finally, by using the symmetric formula, equation ( TTBT) . the matrix 
elements of the higher order terms in the series expansion of the Boltzmann operator are 
calculated to find the matrix elements corresponding to the full time. 

The calculations are done at three different temperatures. Their results are given below. 
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T = 2000°K 
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FIG. 1: 3-D plots and contour plots of the first two terms in the series expansion of the thermal 
propagator at temperature T = 2000 °K for the symmetric Eckart barrier potential. 

1. Symmetric Potential, T = 2000 °K 

In the high temperature limit, the Boltzmann operator is well approximated in terms 
of classical mechanics. Since the TEGA reduces to the classical mechanical Boltzmann 
distribution, one expects it to be accurate in this limit. Figure [1] shows surface and contour 
plots of the matrix elements of the terms K^\f3) and K^ l \f3) at the temperature T = 
2000 °K. In reduced variables, hf3u^ = 1.20 which is small when compared to the reduced 
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T = 2000°K 
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FIG. 2: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator at temperature T = 2000 °K for the symmetric Eckart barrier potential. 

crossover temperature of 2ir between thermal activation and tunneling (w* is the harmonic 
barrier frequency of the Eckart barrier) . From the figure one notes that the first order term 
in the series expansion is indeed very small as compared to the zeroth order term. 

Since it is difficult to quantitatively compare contour plots, in figure EJ a cut of the 
contour plot along the antisymmetric line x' = —x is shown. One notes from the figure, that 
there is virtually no difference between the zero-th and first order results. 



2. Symmetric Potential, T = 200 °K 

Figure [3] shows contour plots of the matrix elements of the imaginary time propagator 
TiZo K U) (£), for iV = 0, 1, 2, 3, respectively at the temperature T = 200 °K (or hfiu x = 12) 
which is below the (reduced) crossover temperature of 2tt. At this temperature, tunneling 
becomes important so that the zero-th order contribution in the series representation is 
no longer sufficient. The contour plot of the matrix elements of the TEGA (zero-th order) 
propagator as shown in panel (a) of the Figure has a single saddle point structure as predicted 
by Liu and Miller 32]. As shown in panel (b) of figure El adding the first order term changes 
the structure of the contour plot completely such that there are two saddle points in the 
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T = 200°K 




5-3-1 1 3 5 



FIG. 3: Contour plots of the matrix elements of the thermal propagator in the TEGA series 
expansion. Panels (a), (b), (c) and (d) correspond to the truncated series expansion of the order 
N = 0, 1, 2, 3 respectively at temperature T = 200 °K for the symmetric Eckart barrier potential. 
In panel (a), the contour values are 10 _n with n = 1, . . . , 7. The highest contour (n = 1) is the solid 
black line, and the lowest contour (n = 7) is the dark blue line. In panel (b) the contour values are 
10 _1 (dashed light blue line), 10" 2 , 10" 3 , 10~ 4 , 10" 5 (solid yellow line), 7.0 x 10~ 6 ,6.0 x 10~ 6 (solid 
red line), 4.0 x 10 -6 (solid light blue line),10 -6 , 10~ 7 (solid dark blue line). The contour values for 
panels (c) and (d) are 10" 1 (dashed red line), 10~ 2 , 10~ 3 , 10~ 4 , 10~ 5 , 6.0 x 10~ 6 (solid yellow line), 
5.0 x 10~ 6 , 4.0 x 10~ 6 (solid red line), 2.5 x 10" 6 (solid light blue line), 10~ 6 , 10~ 7 (solid dark blue 
line) . 
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T = 200°K 
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FIG. 4: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator are shown for the truncated series up to order 3 at temperature T = 200 °K for the 
symmetric Eckart barrier potential. 

graph instead of one. From the figure, it is seen that to obtain the numerically exact results, 
it is necessary to include terms up to order 3. Convergence of the results can also be followed 
from the matrix elements (— x\K N ((3)\x) that are plotted in figure HI 

3. Symmetric Potential, T = 40 °K 

The reduced temperature when T = 40 °K is Hficu* = 60, which is much larger than 2ii, 
so that this temperature corresponds to "deep" tunneling regime. As expected, when the 
temperature is lowered, the calculations become more demanding. The converged results are 
obtained only after including the fifth order term in the series expansion of the imaginary 
time propagator. The convergence can be followed from figure [6] where one dimensional cuts 
along the antisymmetric line for the different terms in the series are plotted. Contour plots 
are given in figure [5] 

Interestingly, the structure becomes even more complicated. Along the antisymmetric 
line, there are three maxima. The saddle point reappears at the origin and there are four 
additional saddle points which are off of the antisymmetry line. This has implications for 
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(a) 



(b) 




FIG. 5: Contour plots of the matrix elements of the zeroth and fifth order truncated series rep- 
resentation of the thermal propagator at the temperature T = 40 °K for the symmetric Eckart 
barrier potential. In panel (a), the contour values are 10 -n with n = 1, ... ,7. The highest con- 
tour (n=l) is the black dashed line, and the lowest contour (n=7) is the dark blue line. In panel 
(6), the contour values are 10 _1 (dashed red line), 10 -2 , 10 -3 , 10 -4 , 10 -5 , 10 -6 (solid yellow line), 
3.0 x 10~ 7 , 10~ 7 , 10~ 8 , -10~ 8 (solid green line), -5.0 x 10~ 8 (solid dark blue line). 



the quantum instanton method 30( , where the two dividing surfaces are taken along the two 
saddle points, as found at the higher temperature. Presumably, one could still take the two 
dividing surfaces to be at the point on the antisymmetric line obtained from the intersection 
of the antisymmetric line and the line that connects each pair of saddle points. However, 
this needs to be studied in more detail. 



4- Asymmetric Potential, T = 200 °K 

The results of the TEGA calculation for the asymmetric potential are shown in figure [7J 
Here too, it was necessary to include all terms up to fifth order for convergence. Due to the 
asymmetry of the barrier, a plot of the matrix elements (— x\K( N > (f3)\x) do not show any 
structure. Besides, it is not helpful for following the convergence of the results, either. For 
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T = 40°K 
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FIG. 6: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator at temperature T = 40 °K for the symmetric Eckart barrier potential. 

these reasons, they are not plotted. In this case, the reduced temperature h{3oj* = 12.5 is 
below the crossover temperature, however not very much lower so that again one has only 
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FIG. 7: Contour plots of the matrix elements of the zeroth and fifth order truncated series repre- 
sentation of the thermal propagator at temperature T = 200 °K for the asymmetric Eckart barrier 
potential. In panel (a), the contour values are 10 _n with n = 0, . . . 6. The highest contour (n = 0) 
appears only on the upper right corner of the graph, and the lowest contour (n = 6) is the green 
line. In panel (b), the contour values are 1, 10" 1 , 10~ 2 , 10~ 3 , 10~ 4 , 3.0 x 10~ 5 , 2.5 x 10~ 5 (yellow 
line), 2.3 x 10" 5 (solid light blue line), 2.1 x 10~ 5 , 10~ 5 (solid dark blue line), 10~ 6 . 

two saddle points. Contour plots of matrix elements for the zeroth order approximation and 
the converged results are given in figure [71 

5. A Comparison Of The Results With The Results Of Miller et. al. 



The matrix elements of the Boltzmann operator for the symmetric and asymmetric Eckart 
barrier had been calculated previously by Miller et. al. However, in that paper, authors 
did not provided contour values in the plots. For this reason, their calculations are dupli- 
cated. First, the Hamiltonian matrix is diagonalized in a sine DVR basis 33]. Then, the 
configuration matrix elements of the Boltzmann operator are calculated by using equation 
1291 In the diagonalization calculation, grid spacing is taken to be 0.05a.u. where the grid 
ranges from — 20a. u. to 20a. u. 
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FIG. 8: A comparison of the results of quantum mechanical and the converged TEGA calculations 
of the matrix elements (—x\K(/3)\x) for the symmetric Eckart barrier potential at temperature 
T = 2000°K . 



For the results of the calculations at temperatures T = 2000 °K and T = 200 °K, when 
the contour plots of the results of the quantum mechanical calculations were prepared (with 
the same contour values that are used in the contour plots of TEGA calculations), there 
was no visual difference between them and the contour plots of the TEGA calculations so 
that they are not given. The results of the quantum mechanical and the TEGA calculations 
of the matrix elements (—%\K(j3)\x) are compared in figures [8] and [9] for the temperatures 
T = 2000 °K and T = 200 °K, respectively. From figure [HI it can be seen that the agreement 
between the quantum mechanical results and the TEGA results is perfect at temperature 
T = 2000 °K. At temperature T = 200 °K, agreement is still quite good, but the results 
differ a little bit in the tunneling region. 

On the other hand, when the results of two different calculations are compared for T = 
40 °K, they differ both quantitively and qualitatively. The contour plot of the results of 
quantum mechanical calculation is given in figure [TUJ In the contour plot, there are still 
two saddle points as in the T = 200 °K case. However, the TEGA calculation predict more 
than two saddle points at that temperature. The difference of the results can also be seen 
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T = 200°K 
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FIG. 9: A comparison of the results of quantum mechanical and the converged TEGA calculations 
of the matrix elements (—x\K(/3)\x) for the symmetric Eckart barrier potential at temperature 
T = 200°K. 

from figure [11] where the results of the two calculations are given for the matrix elements 
(—x\K(J3) \x) . The difference in the results can be attributed to the fact that the way that the 
quantum mechanical calculation is performed is not a proper way of calculating the matrix 
elements of the imaginary time propagator especially at such low temperatures. Because, 
in this calculations one is imposing artificial infinite potential walls at the boundaries which 
does not make sense if the potential surface of the system does not support bound states. 
Since the Eckart potential has a continuous spectrum and do not have any bound states, 
introduction of these artificial infinite walls is a source of error, because it discritizes a 
continuous system. The errors might be expected to be small at high temperatures, which 
can also be seen from the comparisons of the results at temperatures T = 2000 °K and 
T = 200 °K. However, the errors can be significant at low temperatures. Miller et. al. 
did not make these calculations at T = 40 °K anyway. Furthermore, as will be seen in 
section [III Bl increasing number of saddle points is again observed in double well calculations 
both with proper quantum mechanical calculations and also with TEGA and PSTEGA 
calculations. 
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T = 40 °K 




FIG. 10: Contour plot of the results of the quantum mechanical calculation of the matrix elements 
{x2\exp(—/3H)\xi) for the symmetric Eckart barrier potential at temperature T = 40°K. In the 
figure, contour values are 10 -2 (dashed red line), 10 -3 , 10 -4 , 10~ 5 , 10~ 6 , 3 x 10~ 7 (solid light blue 
line), 2 x 10~ 7 , 10~ 7 , 5 x 10~ 8 (solid green line). 



Finally, a quantum mechanical calculation of the matrix elements of the Boltzmann oper- 
ator for the asymmetric Eckart barrier potential at temperature T = 200 °K results in very 
good agreement with the TEGA calculations. The contour plot of the results of quantum 
mechanical calculation is given in figure | 



B. Double Well Potential 



Calculations are done with a one dimensional double well potential surface which has the 
form 

V(x) = V (ax i + bx 2 + c) J (31) 

where Vq = 0.004a.u., a = 1.0(a.u.)~ 4 , b = — 4.0(a.u.)~ 2 , c = 4.0. The mass of the particle is 
taken to be m = 1061. Oa.u.. The calculations are done by using both the bare potential and 
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FIG. 11: A comparison of the results of quantum mechanical and the converged TEGA calculations 
of the matrix elements {—x\K{f3)\x) for the symmetric Eckart barrier potential at temperature 
T = 40°K. 



also the Gaussian averaging. Three different temperatures are used in the calculations which 
are: T = 2000 °K, T = 400 °K, and T = 100 °K. In order to calculate the matrix elements 
of the density matrix, the same numerical procedure that was used in the Eckart Barrier 
calculations is followed. Firstly, an equally spaced grid of 300 points in the range [-3:3] is 
taken. Then, the Coherent States, that are formed around these grid points are propagated 
up to half time using the adaptive step size Cash-Karp Runge-Kutta method 3l|. The 
matrices of the zeroth order approximations to the propagator and the correction operator 
are stored in every time step. Then, the matrices of the higher order terms in the series 
expansion of the propagator is obtained by using the recursion formula (equation E]). Time 
integrations are performed with third order Simpson's integration. Finally, by using the 
symmetric formula, equation fl9]), matrix elements corresponding to full time is calculated. 
Please note that, Gaussian averages of the potential and its derivatives can be calculated 
analytically for the potential surface studied here. Quantum mechanical calculations are 
done with the same parameters given in section IIII A 51 for the Eckart barrier calculations. 
The results that are obtained with the TEGA and the PSTEGA methods will be presented 
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T = 200°K 
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FIG. 12: Contour plot of the results of quantum mechanical calculation of the matrix elements 
{x2\ exp(— f3H)\xx) for the asymmetric Eckart barrier potential at temperature T = 200°K. The 
contour values are the same with the ones that are used in figure [7J 

in the following subsections. Before that, it should be noted that all of the results in this 
section are scaled with the partition function which is Z = Tr[exp(— (3H)]. 

1. T = 2000°K 



In figure [131 a contour plot of the matrix elements of the Boltzmann operator at tem- 
perature T = 2000 °K is shown. That figure is prepared by using the results of the zeroth 
order TEGA calculation with the averaged potential. The results of the other calculations, 
including the quantum mechanical calculation, give almost the identical results so that it 
was not necessary to prepare different graphs for different calculations. This can be realized 
easily from figure [HI in which the matrix elements (— x\ exp(—/3H)\x) are plotted for all of 
the five different calculations. As it can be seen from the figure, there is no way to differen- 
tiate between different graphs. Therefore, at high temperatures use of the bare potential in 
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FIG. 13: A contour plot of the matrix elements of the Boltzmann operator at temperature T = 
2000°K. All of the different calculations gives almost identical results so that the same graph is 
obtained from all of them. The contour values are: 5 x lO^ 1 , 1CT 1 , 5 x 10~ 2 , 10~ 2 , 10~ 3 , 10 -4 , 10~ 5 , 
and 1(T 6 . 

propagation does not lead to a worse zeroth order approximation approximation than the 
use of the Gaussian averaged potential. This can be attributed to the fact that such high 
temperatures involves very short time propagation so that the Gaussians which are very 
narrow initially do not get much broadened. Therefore, use of such narrow Gaussians for 
averaging is like using delta functions, which is the same thing with using bare potentials. 
Consequently, use of the bare potential for such short propagation times does not lead to a 
significant error in the zeroth order approximations. 



2. 



400 °K 



At temperature T = 400 °K, tunneling becomes important so that the zeroth order ap- 
proximations to the thermal propagator does not lead to accurate results. In fact, as it is 
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FIG. 14: One dimensional cuts along the anti-diagonal of the matrix elements of the Boltzmann 
operator is shown at temperature T = 2000°K. In the figure there exists five different plots (results 
of both TEGA and PSTEGA calculations calculated with both the Gaussian averaged potential 
and the bare potential, and also the results of quantum mechanical calculation). However, there 
is no way to differentiate between them. 



shown before 32| they do not even give the correct qualitative picture since they always 
lead to a single saddle point at (0,0). Nevertheless, use of the series representation of the 
propagator converges to the right answers. In order to make a comparison of the TEGA 
and the PSTEGA calculations, the results of the quantum mechanical calculation is shown 
in figure [15l Converged results of the TEGA and PSTEGA calculations are shown in figure 
[TBI For the Gaussian averaged potential calculations , it was necessary to include terms up 
to third order. Convergence of the Gaussian averaged potential calculations are shown for 
the matrix elements (— x\K^ l \f5)\x) in figure [T71 It can be seen from the figure that the 
zeroth order approximation leads to quite bad results both qualitatively and quantitatively. 
Inclusion of the first order correction gives accurate results in some parts of the configuration 
space, but not in the tunneling region. In order to get accurate results also in the tunneling 
region, it is necessary to include even the third order correction terms. 

For the bare potential calculations, it was necessary to include even the forth order 
correction terms in order to converge the results. Besides, it was necessary to use smaller 
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FIG. 15: A contour plot of the matrix elements of the thermal propagator obtained with a quantum 
mechanical calculation at temperature T = 400 °K. The contour values are the same with the ones 
that are used in figure [TUl 

time steps in time integrations. Contour plots of the converged results of the TEGA and 
PSTEGA calculations are shown in figureHHJ In figure [T8| convergence of the results is shown 
for the matrix elements (— x\K^ ((3)\x) for i — 0, . . . , 3. Comparing the results in this figure 
with the results obtained with Gaussian averaged potential calculations, shown in figure [TTJ, 
it can be seen that the use of the bare potential leads to much lower accuracy in the zeroth 
order approximations. Nevertheless, the calculations converge to quite accurate results. 
Besides, despite the fact that fluctuations in the results for the bare potential calculations 
are much bigger than the fluctuations in the results for the Gaussian averaged calculations, 
it can be said that the series representations show similar convergence properties for both of 
the calculations since one of them converges at the third order and the other at the fourth 
order. 
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FIG. 16: Contour plots of the matrix elements of the thermal propagator at temperature T = 
400 °K is shown in panels (a), (b), (c) and (d) for TEGA calculation with the Gaussian averaged 
potential , PSTEGA calculation with the Gaussian averaged potential, TEGA calculation with the 
bare potential and PSTEGA calculation with the bare potential, respectively. The results of the 
Gaussian averaged calculations refer to truncated series of order 3 while the results of the bare 
potential calculations refer to truncated series of order 4. Time step that is used in bare potential 
calculations was half of the time step that is used in Gaussian averaged calculations. The contour 
values are: 5 x 10 _1 (dashed green line), 10" 1 , 10~ 2 , 10~ 3 , 4 x 10~ 4 , 3.2 x 10 -4 (solid yellow line), 
2.7 x 10~ 4 (solid light blue line), 10 -4 , 10 -5 and 10 -6 (solid green line). 
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FIG. 17: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator is shown at temperature T = 400 °K for TEGA and PSTEGA calculations with the 
Gaussian averaged potential. In the figure, panels (a), (6), (c) and (d) refer to the results of the 
truncated series of order N = 0, 1, 2 and 3, respectively. Results of quantum mechanical calculations 
are also shown in each panel. 



3. 



100 °K 



At temperature T = 100°K, tunneling becomes even more important. Besides, the cal- 
culations get more demanding because of the increasing propagation time. 



28 
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FIG. 18: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator at temperature T = 400 °K. The results of TEGA and PSTEGA calculations are 
obtained by using the bare potential. The panels (a), (b), (c) and (d) refer to the results of the 
truncated series of order A^ = 0,1,2, and 3, respectively. The results of quantum mechanical 
calculations are also shown in each panel. 



A contour plot of the results of quantum mechanical calculation is shown in panel (a) of 
figure [T9l Contour plots of the results of TEGA and PSTEGA calculations with the Gaussian 
averaged potentials is shown in panels (b) and (c) of the same figure. At this temperature, 
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FIG. 19: Contour plots of the matrix elements of the Boltzmann operator at temperature 
T = 100 °K. The results of quantum mechanical calculation, TEGA calculation with the Gaussian 
averaged potential, PSTEGA calculation with the Gaussian averaged potential, and TEGA cal- 
culation with the bare potential are shown in panels (a), (b) (c) and (d), respectively. Results of 
the Gaussian averaged calculations refer to truncated series of order 3 while the result of the bare 
potential calculation refers to truncated series of order 4. The contour values are: 5 x 10" 1 (dashed 
red line), 10~\l0~ 2 , 10~ 3 , 5 x 10~ 4 (solid yellow line), 2 x 10~ 4 (solid light blue line), 10~ 4 , 10~ 5 
and 10" 6 . 
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FIG. 20: One dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator at temperature T = 100 °K. The panels (a), (b), (c) and (d) refer to truncated series of 
order N = 1,2,3 and 4, respectively. Results of quantum mechanical calculations are also shown 
in each panel. 

the series for the PSTEGA calculation with the bare potential surface does not converge. 
On the other hand, the series for the TEGA calculation with the bare potential surface still 
converges. A contour plot of the results of the TEGA calculation with the bare potential is 
shown in panel (d) of figure [191 All of the graphs looks very similar. In order to converge the 
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series, it was necessary to include terms up to order 3 for the Gaussian averaged potential 
surface calculations and terms up to order 4 for the bare potential calculations. 

Convergence of the results can again be followed from the antisymmetric line. In figure 
|2"UI one dimensional cuts along the anti-diagonal of the matrix elements of the thermal 
propagator is shown for the truncated series of order N = 1, 2, 3, 4 in panels (a), (b), (c) and 
(d) respectively. The results of quantum mechanical calculation is also shown in each panel. 

Another thing which needs to be noted about the contour plots is the presence of more 
than two saddle points. As in the Eckart Barrier calculations, it is again observed that the 
saddle points move away from the antisymmetric line. From figure UHl it can be seen that 
there exists four saddle points. 

4- A Discussion Of The Results 

In Eckart Barrier calculations, it was observed that the results of the quantum mechanical 
calculations do not agree with the results of TEGA calculations at low temperatures. It 
was argued that the discrepancy between the TEGA and the quantum mechanical results 
should be related with the artificial discretization of a continuous system by imposition of 
wrong boundary conditions to quantum mechanical calculations. On the other hand, double 
well potential surface has a discrete spectrum and do not support any scattering states. 
Therefore, the bound state calculation is a proper way of performing a quantum mechanical 
calculation for calculating the matrix elements of the equilibrium density matrix. The 
agreement of the results of quantum mechanical calculation with the results of the TEGA 
and PSTEGA calculations are very good in this case at all temperatures. This also supports 
that the reason of the discrepancy in the Eckart barrier calculations is related with the 
imposition of wrong boundary conditions to quantum mechanical calculations. 

Considering the zeroth order TEGA and PSTEGA approximations, their accuracy de- 
pends on the temperature. At high temperatures, where the system is almost classical, 
zeroth order approximations lead to accurate results. As the temperature is lowered, accu- 

re quantum effects 



racy of the zeroth order approximations gets worse as expected since t 
becomes important at low temperatures. As shown by Liu and Miller 32|, TEGA always 
leads to a single saddle point at (0,0). A similar analysis can also be made for PSTEGA and 
it can be shown that it is also the case for PSTEGA. Therefore, both TEGA and PSTEGA 
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do not even lead to correct structure at low temperatures where the tunneling effects are 
important. Nevertheless, if the series expansion converges, both of them converge to the 
correct answers even at low temperatures. 

Use of the bare potential results in lower accuracy in the zeroth order approximation com- 
pare to use of the Gaussian averaged potential. Higher accuracy of the results of Gaussian 
averaged calculations can be attributed to the fact that Gaussian averaging of the potential 
surface results from variational principles. Nevertheless, use of the bare potential leads to 
faster integration of equations of motion. However, it also leads to slower convergence of 
the series expansion. Besides, the results fluctuate more during convergence if the bare po- 
tential is used. In this study, it was not possible to converge the results at 100 °K with the 
PSTEGA method if the calculations are done with the bare potential. On the other hand, 
if Gaussian averaged potential is used, the series expansion for the PSTEGA method still 
converges, and it gives accurate results at that temperature. 

One thing needs to be noted about PSTEGA calculations. While integrating the equa- 
tions of motion initial width of the Gaussian wave packet is arbitrary. However, this does not 
mean that one can take any value for the initial width and converge the calculations to the 
correct results. While doing PSTEGA computations, it was necessary to figure out which 
initial width gives the best answers. This is done by comparing the results of PSTEGA 
calculations with the results of TEGA calculations. It was seen that if the initial width of 
the Gaussian is taken to be ~ 1 (in mass weighted coordinates), then the results of TEGA 
and the PSTEGA methods are almost identical for the Gaussian averaged calculations. This 
is true for both the zeroth order approximations and also for the truncated series of any 
order. In other words, n th order PSTEGA expansion and n th order TEGA expansion gives 
identical results within numerical accuracy if the initial arbitrary width of the Gaussians in 
PSTEGA calculations is chosen good. 

IV. DISCUSSIONS AND CONCLUSIONS 

The TEGA and the PSTEGA series representations of the thermal propagator were tested 
for two different potential surfaces. The results show that the number of terms needed in the 
series increases as the temperature is lowered. However, even for a reduced temperature as 
low as hfiuj^ = 60 the expansion converges by the time one reaches the fifth order in the series. 
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In real time, this would make the computation prohibitive, since it would be impossible to 
converge such high order terms using Monte Carlo methods for a multidimensional system. 
In imaginary time, the integrand is much less oscillatory and so there is hope that even when 
dealing with many degrees of freedom, one could converge the higher order terms. 

Even if it turns out that it is not practical to converge the higher order terms of the series 
when the system is "complex" , there is value in the present computation. It does show that 
the series converges rather rapidly and that the series at least in principle does lead to the 
correct result. 

In this paper, numerical convergence properties of the TEGA and the PSTEGA series 
representations of the imaginary time propagator are compared. It is shown that if the initial 
arbitrary width of the Gaussians are chosen good; then, TEGA and PSTEGA methods gives 
identical results within numerical accuracy. Although, the PSTEGA method involves a 
phase space integration, it is shown in the appendix that the momentum coordinates can be 
integrated implicitly, so that the PSTEGA method can also be implemented in configuration 
space. Thus, in both the TEGA and the PSTEGA methods, number of equations of motion 
scales linearly with the dimension of the problem. 

It is seen that the Gaussian averaging is important especially at low temperatures. The 
use of the bare potential leads to very low accuracy for the zeroth order term of the series 
representation such that it causes the series representation not to converge at low tempera- 
tures. 

Another important thing which puts a challenge to semiclassical analysis is that it is 
observed that the number of the saddle points of the matrix elements (x'\K(/3)\x) increases 
as the temperature is lowered. Semiclassically, it is obvious why one should expect two saddle 



points. As analyzed by Miller et. al. 30|, the two saddle points correspond semiclassically 



to the two turning points of the classical periodic orbit on the upside down potential energy 



surface whose half period is h/3 |34j. However, it is found that as the temperature is lowered, 
additional saddle points show up. These point out the need for perhaps a deeper semiclassical 
analysis at low temperature. They also create a challenge to the quantum instanton method 
which used the two saddle points to identify the relevant dividing surfaces for thermal rate 
computations. At the low temperatures, at which one finds more than two saddle points, 
it is not clear which saddle points should be used within the quantum instanton method 
context. This question may become even more acute when dealing with asymmetric systems. 
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An Efficient Way of Integrating Equations of Motion for PSTEGA Calculations 

Equation ( 12TJ) can be integrated implicitly to give 

p(r) = c(r)p(0), (32) 

where c(t) is given by 



c(q, r) = exp I - / dr'rG(r / )" 1 J , (33) 

with the initial condition c(q, 0) = /, which can be integrated with the equation of motion 

^^ = -^G(rrc( q ,r). (34) 

It is useful to define some auxiliary equations of motion that helps to integrate Gaussian 
integrals of p(r). With the following definitions: 

k(r) = f dr'c(q,r') T c(q,T') (35) 

Jo 

w(q,r) = exp (- £ dr'(V(q(r'))) + ^Tr^T 1 ]) (36) 

S(q,r) = l -j T dr'c{q,r')G{r'){VV{q{r'))) (37) 

and integrating the following equations of motion, 

dk(q,r) 

dr 
dw(q,r) 



dr 
ds(q,r) 



c(q,r) T c(q,r), fc(g ,0) = 0, 






(38) 


-w(q,r) ((V(q(T))) + ^Tr[G(rr 1 ]Y 


w(q ,0) = 


= 1, 


(39) 


c(q,T) T G(T)(VV(q(r))), s(q ,0) = 0, 






(40) 



dr 
matrix elements of the zeroth order approximation to the propagator can be obtained as 



(x\K (r)\x') = J ^-(x\k (r)\g(p,q,0))(g(p,q,0)\x') (41) 

dq (2tt)^ 2 / 1 

2n ^m^T) Uet(G(r)G(0)). 



/i;^&T ( * r) (ds^sw))^ 1 ^-^^)- ^ 



where 



l(q, x, x', t) = exp(-i(t(g, r) T k(q, ryH(q, r) + (x - q(r)) T G(r)- 1 (x - q(r)) 

+ ( x '-q )G(0)- 1 (x'-q ))), (43) 

35 



where 



11 1 

%, t) = -s(q, r) + -c(q, r)(x - q(r)) - j{x - q ). (44) 



If the calculations are done with the bare potential, the following term should be added to 
the expression in parenthesis in equations ( )36l) and ( 1391) . 

^Tr[VV T V(q(r))G(r)]; (45) 

and also the Gaussian averagings of the potential and its derivatives are not performed. 

Although, the integration scheme described above increases the number of equations 
of motion per particle. It reduces the phase space integration to a configuration space 
integration so that the total number of equations of motion is greatly reduced. 
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